%% Block FDM
% call this function every timestep
% layer: the current layer No. & medium paras
% model: the current layer's grid
% fields: Vx, Vy, Sxx, Sxy, Syy
function [fields] = calc_blk_field_chrt(layer, model, fields, source)
nx = model.nx;
ny = model.ny;
dx = model.dx;
dy = model.dy;
dt = model.dt;

lambda = layer.lambda;
mu = layer.mu;
rho= layer.rho;
% no_of_lyr=1 means top
% no_of_lyr=n means bottom
no_of_lyr = layer.number;
nlyr = layer.total_lyr;

% the layer where single source is set
src_lyr = source.lyr;
src_vx  = source.vx;
src_vy  = source.vy;
src_sxx = source.sxx;
src_sxy = source.sxy;
src_syy = source.syy;

vx = fields.vx;
vy = fields.vy;
sxx= fields.sxx;
sxy= fields.sxy;
syy= fields.syy;

% Shapiro滤波参数
nu = 0.05; % = 0.20/4
% nu = 0.20/16;

% ghost points
for i = 3:nx-2
    vx(i,2) = 4*vx(i,3) - 6*vx(i,4) + 4*vx(i,5) - vx(i,6);
    vx(i,1) = 4*vx(i,2) - 6*vx(i,3) + 4*vx(i,4) - vx(i,5);
    vy(i,2) = 4*vy(i,3) - 6*vy(i,4) + 4*vy(i,5) - vy(i,6);
    vy(i,1) = 4*vy(i,2) - 6*vy(i,3) + 4*vy(i,4) - vy(i,5);
    vx(i,ny-1) = 4*vx(i,ny-2) - 6*vx(i,ny-3) + 4*vx(i,ny-4) - vx(i,ny-5);
    vx(i,ny) = 4*vx(i,ny-1) - 6*vx(i,ny-2) + 4*vx(i,ny-3) - vx(i,ny-4);
    vy(i,ny-1) = 4*vy(i,ny-2) - 6*vy(i,ny-3) + 4*vy(i,ny-4) - vy(i,ny-5);
    vy(i,ny) = 4*vy(i,ny-1) - 6*vy(i,ny-2) + 4*vy(i,ny-3) - vy(i,ny-4);
end
for j = 3:ny-2
    vx(2,j) = 4*vx(3,j) - 6*vx(4,j) + 4*vx(5,j) - vx(6,j);
    vx(1,j) = 4*vx(2,j) - 6*vx(3,j) + 4*vx(4,j) - vx(5,j);
    vy(2,j) = 4*vy(3,j) - 6*vy(4,j) + 4*vy(5,j) - vy(6,j);
    vy(1,j) = 4*vy(2,j) - 6*vy(3,j) + 4*vy(4,j) - vy(5,j);
    vx(nx-1,j) = 4*vx(nx-2,j) - 6*vx(nx-3,j) + 4*vx(nx-4,j) - vx(nx-5,j);
    vx(nx,j) = 4*vx(nx-1,j) - 6*vx(nx-2,j) + 4*vx(nx-3,j) - vx(nx-4,j);
    vy(nx-1,j) = 4*vy(nx-2,j) - 6*vy(nx-3,j) + 4*vy(nx-4,j) - vy(nx-5,j);
    vy(nx,j) = 4*vy(nx-1,j) - 6*vy(nx-2,j) + 4*vy(nx-3,j) - vy(nx-4,j);
end

% 更新应力分量
for i = 3:nx-2
    for j = 3:ny-2

        dvx_dx = (-vx(i+2,j) + 8*vx(i+1,j) - 8*vx(i-1,j) + vx(i-2,j))/(12*dx);
        dvy_dx = (-vy(i+2,j) + 8*vy(i+1,j) - 8*vy(i-1,j) + vy(i-2,j))/(12*dx);


        dvx_dy = (-vx(i,j+2) + 8*vx(i,j+1) - 8*vx(i,j-1) + vx(i,j-2))/(12*dy);
        dvy_dy = (-vy(i,j+2) + 8*vy(i,j+1) - 8*vy(i,j-1) + vy(i,j-2))/(12*dy);


        sxx(i,j) = sxx(i,j) + dt * ((lambda + 2*mu)*dvx_dx + lambda*dvy_dy);
        syy(i,j) = syy(i,j) + dt * (lambda*dvx_dx + (lambda + 2*mu)*dvy_dy);
        sxy(i,j) = sxy(i,j) + dt * mu * (dvx_dy + dvy_dx);
        % if (check_inside_pml(i,j,nx,ny,no_of_lyr,nlyr) == false)
        %     sxx(i,j) = sxx(i,j) + dt * ((lambda + 2*mu)*dvx_dx + lambda*dvy_dy);
        %     syy(i,j) = syy(i,j) + dt * (lambda*dvx_dx + (lambda + 2*mu)*dvy_dy);
        %     sxy(i,j) = sxy(i,j) + dt * mu * (dvx_dy + dvy_dx);
        % else
        %     pml_para.ui = sxx(i,j);
        %     pml_para.dt = dt;
        %     pml_para.rhs= ((lambda + 2*mu)*dvx_dx + lambda*dvy_dy);
        %     sxx(i,j) = pml_bdry(pml_para,i,j,nx,ny,no_of_lyr,nlyr);
        %     pml_para.ui = sxy(i,j);
        %     pml_para.dt = dt;
        %     pml_para.rhs= (lambda*dvx_dx + (lambda + 2*mu)*dvy_dy);
        %     sxy(i,j) = pml_bdry(pml_para,i,j,nx,ny,no_of_lyr,nlyr);
        %     pml_para.ui = syy(i,j);
        %     pml_para.dt = dt;
        %     pml_para.rhs= mu * (dvx_dy + dvy_dx);
        %     syy(i,j) = pml_bdry(pml_para,i,j,nx,ny,no_of_lyr,nlyr);
        % end

    end
end

% Shapiro滤波
for i = 4:nx-3
    for j = 4:ny-3
        sxx(i,j) = sxx(i,j) + nu * (sxx(i-1,j) - 2*sxx(i,j) + sxx(i+1,j))...
            + nu * (sxx(i,j-1) - 2*sxx(i,j) + sxx(i,j+1));
        sxy(i,j) = sxy(i,j) + nu * (sxy(i-1,j) - 2*sxy(i,j) + sxy(i+1,j))...
            + nu * (sxy(i,j-1) - 2*sxy(i,j) + sxy(i,j+1));
        syy(i,j) = syy(i,j) + nu * (syy(i-1,j) - 2*syy(i,j) + syy(i+1,j))...
            + nu * (syy(i,j-1) - 2*syy(i,j) + syy(i,j+1));
    end
end
% for i = 3:nx-2
%     for j = 3:ny-2
%         sxx(i,j) = sxx(i,j) + ...
%             nu * (sxx(i-2,j) - 4*sxx(i-1,j) + 6*sxx(i,j) - 4*sxx(i+1,j) + sxx(i+2,j))...
%             + nu * (sxx(i,j-2) - 4*sxx(i,j-1) + 6*sxx(i,j) - 4*sxx(i,j+1) + sxx(i,j+2));
%         sxy(i,j) = sxy(i,j) + ...
%             nu * (sxy(i-2,j) - 4*sxy(i-1,j) + 6*sxy(i,j) - 4*sxy(i+1,j) + sxy(i+2,j))...
%             + nu * (sxy(i,j-2) - 4*sxy(i,j-1) + 6*sxy(i,j) - 4*sxy(i,j+1) + sxy(i,j+2));
%         syy(i,j) = syy(i,j) + ...
%             nu * (syy(i-2,j) - 4*syy(i-1,j) + 6*syy(i,j) - 4*syy(i+1,j) + syy(i+2,j))...
%             + nu * (syy(i,j-2) - 4*syy(i,j-1) + 6*syy(i,j) - 4*syy(i,j+1) + syy(i,j+2));
%     end
% end

% set the upper interface of 1st layer as free surf.
if no_of_lyr == 1
    top_free_surf_j = 3;
    syy(:,top_free_surf_j) = 0;
    sxy(:,top_free_surf_j) = 0;
end

% add the source all around the grid block
if no_of_lyr == src_lyr
    vx = vx + src_vx;
    vy = vy + src_vy;
    sxx= sxx+ src_sxx;
    sxy= sxy+ src_sxy;
    syy= syy+ src_syy;
end

% ghost points
for i = 3:nx-2
    sxx(i,2) = 4*sxx(i,3) - 6*sxx(i,4) + 4*sxx(i,5) - sxx(i,6);
    sxx(i,1) = 4*sxx(i,2) - 6*sxx(i,3) + 4*sxx(i,4) - sxx(i,5);
    sxy(i,2) = 4*sxy(i,3) - 6*sxy(i,4) + 4*sxy(i,5) - sxy(i,6);
    sxy(i,1) = 4*sxy(i,2) - 6*sxy(i,3) + 4*sxy(i,4) - sxy(i,5);
    syy(i,2) = 4*syy(i,3) - 6*syy(i,4) + 4*syy(i,5) - syy(i,6);
    syy(i,1) = 4*syy(i,2) - 6*syy(i,3) + 4*syy(i,4) - syy(i,5);
    sxx(i,ny-1) = 4*sxx(i,ny-2) - 6*sxx(i,ny-3) + 4*sxx(i,ny-4) - sxx(i,ny-5);
    sxx(i,ny) = 4*sxx(i,ny-1) - 6*sxx(i,ny-2) + 4*sxx(i,ny-3) - sxx(i,ny-4);
    sxy(i,ny-1) = 4*sxy(i,ny-2) - 6*sxy(i,ny-3) + 4*sxy(i,ny-4) - sxy(i,ny-5);
    sxy(i,ny) = 4*sxy(i,ny-1) - 6*sxy(i,ny-2) + 4*sxy(i,ny-3) - sxy(i,ny-4);
    syy(i,ny-1) = 4*syy(i,ny-2) - 6*syy(i,ny-3) + 4*syy(i,ny-4) - syy(i,ny-5);
    syy(i,ny) = 4*syy(i,ny-1) - 6*syy(i,ny-2) + 4*syy(i,ny-3) - syy(i,ny-4);
end
for j = 3:ny-2
    sxx(2,j) = 4*sxx(3,j) - 6*sxx(4,j) + 4*sxx(5,j) - sxx(6,j);
    sxx(1,j) = 4*sxx(2,j) - 6*sxx(3,j) + 4*sxx(4,j) - sxx(5,j);
    sxy(2,j) = 4*sxy(3,j) - 6*sxy(4,j) + 4*sxy(5,j) - sxy(6,j);
    sxy(1,j) = 4*sxy(2,j) - 6*sxy(3,j) + 4*sxy(4,j) - sxy(5,j);
    syy(2,j) = 4*syy(3,j) - 6*syy(4,j) + 4*syy(5,j) - syy(6,j);
    syy(1,j) = 4*syy(2,j) - 6*syy(3,j) + 4*syy(4,j) - syy(5,j);
    sxx(nx-1,j) = 4*sxx(nx-2,j) - 6*sxx(nx-3,j) + 4*sxx(nx-4,j) - sxx(nx-5,j);
    sxx(nx,j) = 4*sxx(nx-1,j) - 6*sxx(nx-2,j) + 4*sxx(nx-3,j) - sxx(nx-4,j);
    sxy(nx-1,j) = 4*sxy(nx-2,j) - 6*sxy(nx-3,j) + 4*sxy(nx-4,j) - sxy(nx-5,j);
    sxy(nx,j) = 4*sxy(nx-1,j) - 6*sxy(nx-2,j) + 4*sxy(nx-3,j) - sxy(nx-4,j);
    syy(nx-1,j) = 4*syy(nx-2,j) - 6*syy(nx-3,j) + 4*syy(nx-4,j) - syy(nx-5,j);
    syy(nx,j) = 4*syy(nx-1,j) - 6*syy(nx-2,j) + 4*syy(nx-3,j) - syy(nx-4,j);
end

% 更新速度分量
for i = 3:nx-2
    for j = 3:ny-2

        dsxx_dx = (-sxx(i+2,j) + 8*sxx(i+1,j) - 8*sxx(i-1,j) + sxx(i-2,j))/(12*dx);
        dsxy_dx = (-sxy(i+2,j) + 8*sxy(i+1,j) - 8*sxy(i-1,j) + sxy(i-2,j))/(12*dx);

        dsxy_dy = (-sxy(i,j+2) + 8*sxy(i,j+1) - 8*sxy(i,j-1) + sxy(i,j-2))/(12*dy);
        dsyy_dy = (-syy(i,j+2) + 8*syy(i,j+1) - 8*syy(i,j-1) + syy(i,j-2))/(12*dy);


        vx(i,j) = vx(i,j) + dt/rho * (dsxx_dx + dsxy_dy);
        vy(i,j) = vy(i,j) + dt/rho * (dsxy_dx + dsyy_dy);
        % if (check_inside_pml(i,j,nx,ny,no_of_lyr,nlyr) == false)
        %     vx(i,j) = vx(i,j) + dt/rho * (dsxx_dx + dsxy_dy);
        %     vy(i,j) = vy(i,j) + dt/rho * (dsxy_dx + dsyy_dy);
        % else
        %     pml_para.ui = vx(i,j);
        %     pml_para.dt = dt;
        %     pml_para.rhs= (dsxx_dx + dsxy_dy)/rho; % right hand side 方程右端项
        %     vx(i,j) = pml_bdry(pml_para,i,j,nx,ny,no_of_lyr,nlyr);
        %     pml_para.ui = vy(i,j);
        %     pml_para.dt = dt;
        %     pml_para.rhs= (dsxy_dx + dsyy_dy)/rho; % right hand side 方程右端项
        %     vy(i,j) = pml_bdry(pml_para,i,j,nx,ny,no_of_lyr,nlyr);
        % end
    end
end

% Shapiro滤波
for i = 4:nx-3
    for j = 4:ny-3
        vx(i,j) = vx(i,j) + nu * (vx(i-1,j) - 2*vx(i,j) + vx(i+1,j))...
            + nu * (vx(i,j-1) - 2*vx(i,j) + vx(i,j+1));
        vy(i,j) = vy(i,j) + nu * (vy(i-1,j) - 2*vy(i,j) + vy(i+1,j))...
            + nu * (vy(i,j-1) - 2*vy(i,j) + vy(i,j+1));
    end
end
% for i = 3:nx-2
%     for j = 3:ny-2
%         vx(i,j) = vx(i,j) + ...
%             nu * (vx(i-2,j) - 4*vx(i-1,j) + 6*vx(i,j) - 4*vx(i+1,j) + vx(i+2,j))...
%             + nu * (vx(i,j-2) - 4*vx(i,j-1) + 6*vx(i,j) - 4*vx(i,j+1) + vx(i,j+2));
%         vy(i,j) = vy(i,j) + ...
%             nu * (vy(i-2,j) - 4*vy(i-1,j) + 6*vy(i,j) - 4*vy(i+1,j) + vy(i+2,j))...
%             + nu * (vy(i,j-2) - 4*vy(i,j-1) + 6*vy(i,j) - 4*vy(i,j+1) + vy(i,j+2));
%     end
% end

fields.vx = vx;
fields.vy = vy;
fields.sxx= sxx;
fields.sxy= sxy;
fields.syy= syy;

end